program Main

    use omp_lib
    use Globals
    use ModuleSolveEqm
    use ModuleQERRORS
    use ModuleElasticities

    implicit none
    
  
    ! Local variable declarations
    real(8) :: param_vec_init(3)    ! vector of tax function parameters (initial system): input for function solving for equilibrium given taxes
    real(8) :: param_vec_lr(3)      ! vector of tax function parameters (final system): input for function solving for equilibrium given taxes
    real(8) :: welfare              ! auxiliary variable to store welfare numbers
    integer :: clck_counts_beg, clck_counts_end, clck_rate  ! for timing (true time; not CPU time with parallelization)
    real(8) :: ParamCalib(5)        ! calibrated parameters
    real(8) :: eqm_results(5)       ! equilibrium values from calibration
    real(8) :: V_TR_init(NS)        ! value function in first period of transition
    
    ! Start timer
    call system_clock ( clck_counts_beg, clck_rate )
    
    ! Choose normalization version
    norm_version = 2            ! indicator for choice of normalization: 1) mean income; 2) median income; 3) mean income of the third quintile
     
    ! Calibrated parameters
    open(11,  file = 'Input/ParamCalib.txt',      status = 'unknown')
    read(11, *) ParamCalib
    close(11)
    beta        = ParamCalib(1)             ! discount factor
    B           = ParamCalib(2)             ! labor disutility parameter
    G           = ParamCalib(3)             ! government spending
    Debt        = ParamCalib(4)             ! government debt
    ymean       = ParamCalib(5)             ! mean income
    ymeanCalib  = ParamCalib(5)             ! mean income
    
    ! Set tax function parameters
    open(11,  file = 'Input/eqm_results.txt',      status = 'unknown')
    read(11, *) eqm_results
    close(11)
    param_vec_init(1) = eqm_results(3)     ! gamma
    param_vec_init(2) = eqm_results(2)     ! lambda
    param_vec_init(3) = eqm_results(5)     ! rt
    omegat = 0d0                           ! omegat
    print *, 'Tax parameters initially: gamma, lambda, rt = ', param_vec_init
    print *, ''
    
    ! Guesses for equilibrium objects
    mt     = eqm_results(4)             ! mt
    r      = eqm_results(1)             ! interest rate
    lambda = param_vec_init(2)          ! lambda
    print*, 'Guesses for r = ', r,' and lambda = ', lambda
    print *, ''

    ! Solve for initial steady state
    print *, 'Computing initial steady-state'
    print *, ''
    
    ! Solve for equilibrium given tax function parameters
    welfare = welfare_eqm(param_vec_init,1)

    ! Store eqm objects
    r_init = r                  ! initial equilibrium interest rate
    mt_init = mt                ! initial equilibrium level transfer
    error_vec_init = error_vec  ! initial equilibrium vector of errors
    Kagg_init = Kagg            ! initial equilibrium capital stock

    ! Solve for new steady state
    print *, ''
    print *, 'Computing new steady-state'
    print *, ''

    ! Set tax function parameters for final equilibrium
    param_vec_lr(1) = eqm_results(3)  ! gamma
    param_vec_lr(2) = eqm_results(2)  ! lambda
    param_vec_lr(3) = eqm_results(5)  ! rt

    ! Solve for equilibrium given tax function parameters
    welfare = welfare_eqm(param_vec_lr,2)

    ! Store eqm objects
    r_lr = r                    ! new equilibrium interest rate
    lambda_lr = lambda          ! new equilibrium tax function level
    mt_lr = mt                  ! new equilibrium level transfer
    error_vec_lr  = error_vec   ! new equilibrium vector of errors

    ! Solving for transition
    print *, ''
    print *, 'Start transition'
    print *, ''

    
    ! Guesses for transition
    r_TR = r_lr                                     ! interest rate path
    lambda_TR = lambda_lr                           ! tax level path
    mt_TR = mt_lr                                   ! transfer level path
    wge_TR = (1.0D0-alpha)/(r_TR+delta)             ! wage path
    wge_TR = alpha*(wge_TR**((1d0-alpha)/alpha))      ! wage path

    epsilon = 0.0D0
    epsilon(1) = 0.01d0 ! Small change in tax, period 1 only

    call Model_VF_TR


    call elasticities

    ! End timer
    call system_clock ( clck_counts_end, clck_rate )
    print *, 'Time in seconds = ' , (clck_counts_end - clck_counts_beg) / real (clck_rate)

end program Main
